Merger times of double compact object systems

The timescale for the merger due to gravitational wave emission of a double compact object (DCO) binary can be analytically approximated by

\begin{aligned} t_\mathrm{merger} & = \dfrac{15}{304} \dfrac{a_{0}^{4} c^{5}}{G^{3} m_{1} m_{2} (m_{1} + m_{2})} \times \left[\left(1+e_{0}^{2}\right) e_{0}^{-12/19} \left(1 + \dfrac{121}{304} e_{0}^{2} \right)^{-870/2299}\right]^{4} \times \int_{0}^{e_{0}} f(e) \, \mathrm{d}e \end{aligned}

where $c$ is the speed of light, $G$ the gravitational constant, $m_1$ and $m_2$ are the masses of the compact objects and $a_0$ and $e_0$ are the separation and eccentricity of the DCO binary at the formation of the system, respectively. Additionally,

\begin{aligned} f(e) & = \dfrac{e^{29/19} \left[1+\left(121/304\right)e^{2}\right]^{1181/2299}}{\left(1-e^{2}\right)^{3/2}} \end{aligned}

such that the integration proceeds from the formation of the DCO binary until its merger. For a detailed derivation, see Peters 1964.

From the expression of the $t_\mathrm{merger}$, we see that given a set of masses ($m_1$, $m_2$) there is limit in the binary parameters ($a$, $e$) such that they merge in less than the Hubble time.

In order to find this region, we define a function f as

\begin{equation} f = t_\mathrm{merger}(m_1, m_2, a, e) - t_\mathrm{Hubble} \end{equation}

Then, given a set of masses $m_1$ and $m_2$, one can find the values of $a_0$ and $e_0$ which are roots of the function $f$,

\begin{equation} f = t_\mathrm{merger}(m_1, m_2, a_0, e_0) - t_\mathrm{Hubble} = 0 \end{equation}

With a call to a root solver,like the bisect method of scipy.optimize, we can find this region, which we show below.